/* -*- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -*- */
// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
/* tests/test-integer-tools.C
 * Copyright (C) 2011 LinBox
 *
 * Written by Brice Boyer <Brice.Boyer@imag.fr>
 *
 * See COPYING for license information.
 */


/*! @file  tests/test-integer-tools.C
 * @ingroup tests
 * Tests  integer tools such as Matrix-Matrix multiplication, Matrix-vector, p-adic
 * conversions,...
 * @test Tests  integer tools such as Matrix-Matrix multiplication, Matrix-vector, p-adic
 * conversions,...
 */

#include <iostream>
#include <cassert>
#include "linbox/linbox-config.h"
#include "test-common.h"
#include "linbox/util/timer.h"

#include "linbox/algorithms/integer-tools.h"
#include "linbox/matrix/random-matrix.h"

using namespace LinBox ;

/*! Check integer matrix multiplication D = alpha A B + beta C.
 * We only perform dot-products on random rows/columns
 */
template<class Matrix, class Field>
bool FastCheckMatMul(const Field                   & F,
		     const Matrix                  & D,
		     const typename Field::Element & alpha,
		     const Matrix                  & A,
		     const Matrix                  & B,
		     const typename Field::Element & beta,
		     const Matrix                  & C
		 )
{
	size_t M = C.rowdim();
	size_t N = C.coldim();
	size_t K = A.coldim();
	if(!A.rowdim() || !A.coldim() || !B.coldim())
		return (!C.rowdim()||!C.coldim());

	typedef typename Field::Element Element ;

	for (size_t l = 0 ; l < (size_t) std::log(M*N) ; ++l) {
		size_t i = Integer::random((Integer)M);
		size_t j = Integer::random((Integer)N);

#ifdef _LB_DEBUG
		// std::cout << "checking " << i << ',' << j << std::endl;
#endif


		Element Cij ;
		F.init(Cij,0UL);
		for (size_t k = 0 ; k < K ; ++k)
			F.axpyin(Cij,A.getEntry(i,k),B.getEntry(k,j));
		if (!(F.areEqual(alpha*Cij+beta*C.getEntry(i,j),D.getEntry(i,j)))) {
			// #ifdef _LB_DEBUG
			std::cout << "@(" <<i << ',' << j << ") : " <<  alpha << '*' << Cij << "+" << beta << '*' << C.getEntry(i,j) << "!=" << D.getEntry(i,j) << std::endl;
			// #endif
			return false ;
		}
	}

	return true;

}

bool RandomCheckMatMul(
		       const BlasMatrix<Integer>   & DD
		       , const Integer             & alpha
		       , const BlasMatrix<Integer> & A
		       , const BlasMatrix<Integer> & B
		       , const Integer             & beta
		       , const BlasMatrix<Integer> & C
		 )
{
	PID_integer ZZ ;
	BlasMatrix<Integer> D(DD.rowdim(),DD.coldim());
	MatrixDomain<PID_integer> ZMD(ZZ);
	ZMD.muladd(D,beta,C,alpha,A,B);


	size_t M = C.rowdim();
	size_t N = C.coldim();
	// size_t K = A.coldim();
	if(!A.rowdim() || !A.coldim() || !B.coldim())
		return (!C.rowdim()||!C.coldim());

	typedef Modular<double>       ModularField ;
	typedef ModularField::Element      Element ;

	Integer p = Integer::random_between(10,12) ;
	ModularField Zp (p);

	BlasMatrixDomain<ModularField> BMD (Zp);
	MatrixDomain<ModularField>      MD (Zp);

	BlasMatrix<Element> Ap(A,Zp);
	BlasMatrix<Element> Bp(B,Zp);
	BlasMatrix<Element> Cp(C,Zp);
	BlasMatrix<Element> Dp(D,Zp);
	BlasMatrix<Element> Ep(M,N);

	Element ap, bp ;
	Zp.init(ap,alpha);
	Zp.init(bp,beta);

	BMD.muladd(Ep,bp,Cp,ap,Ap,Bp);

	bool pass = MD.areEqual(Ep,Dp);
	if (!pass) {
		// if (MD.areEqual(E,D))
			// std::cout << " *** ERROR in muladd *** " << std::endl;
		// else
			// E.write(std::cout << "EEEE :=", true) << ";" << std::endl;
		std::cout << "#########################################" << std::endl;
		std::cout << "p := " << p << ';' << std::endl;
		std::cout << "ap,bp := " << ap << ',' << bp << ';' << std::endl;
		Ap.write(std::cout << "Ap :=", true) << ";" << std::endl;
		Bp.write(std::cout << "Bp :=", true) << ";" << std::endl;
		Cp.write(std::cout << "Cp :=", true) << ";" << std::endl;
		Dp.write(std::cout << "Dp :=", true) << ";" << std::endl;
		Ep.write(std::cout << "Ep :=", true) << "; (bad)" << std::endl;
		std::cout << "alpha,beta := " << alpha << ',' << beta << ';' << std::endl;
		A.write(std::cout << "A :=",true) << ';' << std::endl;
		B.write(std::cout << "B :=",true) << ';' << std::endl;
		C.write(std::cout << "C :=",true) << ';' << std::endl;
		D.write(std::cout << "E :=",true) << ';' << std::endl;
		std::cout << "evalm(E-alpha*A.B-beta*C);" << std::endl;
		std::cout << "#########################################" << std::endl;

	}

	MD.muladd(Ep,bp,Cp,ap,Ap,Bp);
	pass = MD.areEqual(Ep,Dp);
	if (!pass) {
		std::cout << " *** ERROR *** " << std::endl;
	}

	return pass ;

}

bool test_matmul(int m, int n, int k,
		 const Integer & alpha,
		 const Integer & beta)
{
	bool pass = true ;
	typedef BlasMatrix<Integer> IntMat ;
	IntMat A(m,k);
	IntMat B(k,n);
	IntMat C(m,n);
	IntMat D(m,n);

	typedef RandomIntegerIter<>                       IntRandIter ;
	typedef RandomDenseMatrix<IntRandIter, PID_integer > IntRand_t;

	size_t lA = 22 ;
	size_t lB = 23 ;
	size_t lC = 24 ;

	IntRandIter RA(lA);
	IntRandIter RB(lB);
	IntRandIter RC(lC);

	PID_integer ZZ ;

	IntRand_t Arand (ZZ,RA);
	IntRand_t Brand (ZZ,RB);
	IntRand_t Crand (ZZ,RC);


	Arand.random<IntMat>(A);
	Brand.random<IntMat>(B);
	Crand.random<IntMat>(C);


	// D = alpha A B + beta C ;
	MMmul_naif(D,alpha,A,B,beta,C);
	pass = RandomCheckMatMul(D,alpha,A,B,beta,C);
	if (!pass) {
		std::cout << "naive mul failed" << std::endl;
		std::cout << "alpha,beta := " << alpha << ',' << beta << ';' << std::endl;
		A.write(std::cout << "A :=",true) << ';' << std::endl;
		B.write(std::cout << "B :=",true) << ';' << std::endl;
		C.write(std::cout << "C :=",true) << ';' << std::endl;
		D.write(std::cout << "E :=",true) << ';' << std::endl;
		std::cout << "evalm(E-alpha*A.B-beta*C);" << std::endl;

	}
	MMmul_naif_in(C,alpha,A,B,beta);
	MatrixDomain<PID_integer> MD(ZZ);
	bool pass2 = MD.areEqual(C,D);
	if (!pass2) {
		std::cout << "naive mulin failed " << std::endl;
		pass |= pass2;
		std::cout << "alpha,beta := " << alpha << ',' << beta << ';' << std::endl;
		A.write(std::cout << "A :=",true) << ';' << std::endl;
		B.write(std::cout << "B :=",true) << ';' << std::endl;
		C.write(std::cout << "C :=",true) << ';' << std::endl;
		D.write(std::cout << "E :=",true) << ';' << std::endl;
		std::cout << "evalm(E-C);" << std::endl;
	}
	return !pass ;

}

int main(int ac, char ** av)
{
#if 0
	/*  Argument parsing/setting */

	static size_t       n = 50;    /*  Taille */
	static size_t       m = 50;    /*  Taille */
	static size_t       p = 22;    /*  PrimeSize */
	// static size_t    seed =  0;    /*  ! unused */
	static size_t   iters = 20;    /* _LB_REPEAT */

        static Argument as[] = {
                { 'n', "-n N", "Set number of primes.", TYPE_INT , &n },
                { 'p', "-p P", "Set size of test primes.", TYPE_INT , &p },
                { 'i', "-i I", "Perform each test for I iterations.",     TYPE_INT, &iters },
		END_OF_ARGUMENTS
        };

	parseArguments (ac, av, as);
#endif

	bool pass = true ;

	srand(time(NULL));             // seeding

	size_t m,n,k ;
	m = 10 ; n = 10 ; k = 10 ;
	Integer alpha, beta ;
	alpha = 0 ; beta = 1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = 0 ; beta = Integer::random() ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = 1 ; beta = 0 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = 1 ; beta = 1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = 1 ; beta = -1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = 1 ; beta = Integer::random() ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = -1 ; beta = 0 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = -1 ; beta = 1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = -1 ; beta = -1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = -1 ; beta = Integer::random() ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = Integer::random() ; beta = 0 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha =  Integer::random() ; beta = 1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = Integer::random() ; beta = -1 ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);
	alpha = Integer::random() ; beta = Integer::random() ;
	std::cout << alpha << ',' << beta << std::endl;
	pass |= test_matmul(m,n,k,alpha,beta);

	return !pass ;


}
